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Afcifract-Optimization  of  similarity  metrics  is  an  important 
component  of  many  biomedical  image  registration  applications 
that  greatly  affects  registration  outcome.  Selection  of  an 
optimization  approach  is  dependent  on  the  modalities  being 
registered,  the  characteristics  of  the  similarity  function,  and 
available  computational  resources.  This  paper  addresses 
optimization  approaches  for  2D  and  3D  rigid  body  intra-  and 
multimodal  ultrasound  registration.  Stochastic  and  direct 
techniques  are  compared  for  mutual  information  and 
correlation  ratio  functions.  A  new  direct/stochastic  hybrid 
approach  based  on  the  tabu  search  is  also  proposed. 
Visualization  and  experimental  results  suggest  the  usefulness 
of  such  an  approach.  Results  also  show  that  the  tabu  hybrid 
technique  compares  favorably  with  traditional  techniques, 
based  on  accuracy  and  number  of  function  evaluations. 
Keywords  -  Image  registration,  ultrasound,  search, 
optimization 

I.  Introduction 

Rigid  body  registration  of  2D-to-2D  MRI  and  CT 
biomedical  images  can  be  performed  efficiently  and 
accurately  with  well-established  methods.  However,  2D-to- 
3D  and  3D-to-3D  ultrasound  (US)  registration  and 
subsequent  fusion  is  a  relatively  new  and  challenging  field, 
as  speckle  and  low  signal -to-noise  ratio  (SNR)  make  both 
human  and  computerized  analysis  difficult.  There  are 
generally  three  major  components  in  image  registration:  (1) 
definition  of  a  search  space;  (2)  selection  of  a  similarity 
measure;  and  (3)  choice  of  an  optimization  approach  [1,2]. 
Rigid  body,  affine,  and  elastic  transformations  are  common 
in  US  registration.  Mutual  information  (MI)  and  its  variants 
[3-5]  and  the  correlation  ratio  (CR)  [6,7]  have  been 
proposed  as  robust  similarity  metrics  for  biomedical  images, 
including  US  B-scans.  This  paper  compares  optimization 
methods  for  US  registration  using  these  two  functions.  The 
global  maximum  of  similarity  functions  is  generally 
considered  to  correspond  to  correct  registration.  This 
assumption  is  made  in  the  current  study,  although  it  is  not 
true  in  all  cases,  and  some  degree  of  human  intervention  is 
often  required. 

An  adaptation  of  the  continuous  tabu  search,  combined 
with  a  direct  approach,  is  presented  as  an  accurate  hybrid 
method  for  global  optimization.  The  search  domain  is  the 
space  of  rigid  body  transformations. 

II.  Stochastic  and  Direct  Search  Techniques 

Many  direct  (deterministic)  approaches  are  currently 
used  to  optimize  registration  functionals,  including  those 
requiring  only  function  evaluations  (e.g.  Nelder-Mead 


simplex)  and  those  requiring  derivative  information,  such  as 
gradient  descent,  and  Powell-type  algorithms  [2, 3, 5-7]. 
They  are  primarily  local  optimizers  that  often  rapidly 
converge  to  an  exact  local  optimum.  However,  especially 
for  multiextremal  functions,  they  tend  to  converge  to  the 
first  identified  local  optima  [8],  Continuity,  differentiability, 
and  convexity  constraints  are  placed  on  the  functionals. 
Stochastic  global  optimizers  and  stochastic-direct  hybrids 
are  often  used  to  address  these  problems.  These  methods, 
which  include  genetic  algorithms,  simulated  annealing,  and 
clustering  methods  [9],  can  converge  to  global  optima  of 
complicated  functions.  They  are  robust,  can  quickly  locate 
global  optima,  and  usually  do  not  “get  trapped"  in  local 
extrema.  However,  because  they  are  probabilistic,  they 
sometimes  do  not  locate  the  true  global  optimum,  but  only 
nearby  solutions.  The  search  may  be  time-expensive 
because  of  the  combinatorial  nature  of  sampling 
multidimensional  space.  Nonetheless,  stochastic  methods 
have  proven  effective  for  biomedical  image  registration, 
where  genetic  algorithms  [10,11]  and  simulated  annealing 
[10,12]  have  been  used.  Another  global  optimizer,  the  tabu 
search  (TS),  has  been  successfully  used  to  solve  many 
discrete  optimization  problems  [13],  The  current  work 
presents  an  adaptation  of  the  continuous  tabu  search  in 
which  a  likely  area  for  the  global  optimum  is  determined, 
and  Powell’s  direction  set  algorithm  (PDS)  with  Brent’s  line 
optimization  [14]  identifies  the  global  optimum.  Continuous 
tabu  search  techniques  have  recently  appeared  in  the 
literature  [15-17],  but,  thus  far,  applications  to  image 
registration  have  not  been  reported.  Results  are  compared 
with  those  of  another  powerful  technique,  adaptive 
simulated  annealing  (ASA)  [18],  also  using  PDS  to  find  the 
global  optimum. 

III.  Hybrid  Continuous  Tabu  Search 

The  TS  is  an  iterative  combinatorial  technique. 
Application-specific  heuristics  are  implemented  on  top  of  a 
basic  structure.  One  underlying  concept  of  any  TS  approach 
is  that  “cycles"  are  avoided  while  searching  for  a  potential 
new  optimum  value  by  placing  previously  visited  solutions 
into  a  finite  length  first-in-first-out  “tabu  list”.  The  current 
study  adapts  the  approach  in  [16]  to  image  registration.  The 
proposed  tabu  hybrid  (TH)  consists  of  the  following  stages; 

1.  The  search  space  is  discretized  [15].  Discretization 
provides  effective  neighborhood  generation  and 
facilitates  tabu  list  and  promising  list  searches  (see 
below).  In  the  current  study,  each  parameter  is 
partitioned  into  100  points. 
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2.  In  the  diversification  stage,  k  neighbors  of  the  current  n- 
D  point  s  are  generated  at  each  iteration  by  randomly 
selecting  j  <  n  dimensions  of  s,  and  perturbing  the 
variables  in  those  dimensions.  A  neighborhood  B( s;  r) 
is  a  set  of  nested  “balls”  centered  on  the  current 
solution  s,  with  maximum  radius  r.  It  contains  all 
points  s'  such  that  ||s'  -  sA||  <  r,  where  ||- 1|  denotes  a 
Euclidean  norm.  The  set  of  balls  centered  on  s  with 
radius  h0,  hi,  ...,  hk  is  defined  as  [16]: 

B,( s;  Ki,  ^ )  =  [s'  I  Ki  <  ||s'  -  sA]|  <  hi},  i  =1  (1) 

3.  The  neighbor  with  the  best  function  value  is  selected  as 
the  new  starting  point  if  it  is  not  within  a  neighborhood 
of  a  point  on  the  tabu  list.  If  all  k  neighbors  degrade  the 
functional  more  than  a  threshold,  the  current  point  is 
added  to  a  promising  list.  The  current  point  is  always 
added  to  the  tabu  list.  The  process  continues  until  a 
specified  number  of  iterations  without  a  new  promising 
area  are  reached  (20  in  this  study),  or  until  a  maximum 
number  of  iterations  is  performed  (here,  100). 

4.  Also  in  the  diversification  stage,  the  search  is  restarted 
from  a  randomly  chosen  reflection  of  the  original 
(initial)  solution  if  function  evaluations  do  not  improve 
beyond  the  best  value  obtained,  after  a  fixed  number  of 
iterations  (20  in  this  study). 

5.  An  “affine  shaker”  algorithm  [15]  is  applied  to  locate 
the  best  point  in  each  promising  area.  The  shape  of  the 
promising  area  is  adapted  during  the  search  to  include 
areas  of  the  search  space  that  may  have  been  missed 
during  diversification. 

6.  A  direct  local  optimizer  performs  an  intensified  search 
on  the  best  point  from  step  (5).  Powell’s  algorithm  was 
chosen  because  of  its  accuracy  and  speed. 

The  tabu  search  is  not  a  “brute  force”  approach.  At  each 
iteration,  the  point  with  the  lowest  function  value  in  the 
neighborhood  is  selected  as  the  new  point;  it  need  not  be 
better  than  the  current  point.  In  this  way,  promising  regions 
where  the  global  optimum  may  exist  are  identified,  and 
local  optima  are  escaped.  The  tabu  list  prevents  cyclic 
searching,  and  points  on  the  promising  list  are  explored 
more  intensively  with  the  affine  shaker  and  PDS.  In  some 
implementations,  a  promising  list  is  not  used  [17],  Search 
parameters  can  also  be  tuned  to  increase  performance. 


IV.  Registration  Metrics 
Normalized  mutual  information  is  given  as  [4]: 


I(At,Bt ) 


H(At)  +  H(Bt ) 
H(At,Bt ) 


(2) 


where  H(AT)  denotes  the  entropy  of  estimate  of  image  A 
under  transformation  T,  and  H(Ar,  Br)  denotes  the  joint 
entropy  estimate  of  images  AT  and  BT.  MI  is  commutative. 


(a)  (b) 

Fig.  1 .  Similarity  metrics  for  (a)  MI  and  (b)  CR  for  registration  with  a  high 
SNR  shoulder  B-scan. 


y-trans. 


x-trans. 


This  metric  is  invariant  under  the  size  of  image-to-volume 
overlap  [4],  and  assumes  values  from  0  to  2. 

The  correlation  ratio,  r\(A7]BT),  is  given  as  [7]: 


_  Var[E( Ar  I BT )] 
Var(A, ) 


(3) 


This  metric,  assuming  values  from  0  to  1,  can  also  be 
generalized  for  gray  level  mapping  and  for  weighting 
interpolation  values  [6],  In  general,  A  r|(/i7l,47). 

In  Fig.  1,  plots  of  MI  and  CR  as  functions  of  translation 
are  shown  for  a  high  SNR  B-scan  of  a  shoulder.  The  strong 
differences  in  the  functions  suggest  that  the  metrics  may 
produce  different  registration  results,  and  should  be 
considered  in  selecting  an  optimization  method. 


V.  Experiments 

To  empirically  assess  the  differences  in  the  functionals 
for  B-scan  registration,  MI  and  CR  surfaces  were  generated 
as  functions  of  translation  and  rotation  (see  [3]  and  [7]  for  a 
theoretical  discussion).  50x50-pixel  regions  of  reference 
images  were  registered  with  the  reference.  This  idealized 
situation  provides  an  indication  of  the  characteristics  of  US 
similarity  metrics.  Functionals  were  generated  for:  (1)  A 
region  of  a  low  SNR  kidney  phantom  (simulated  with  Field 
II  Matlab  software  [19])  registered  to  the  original  896x512 
image,  shown  in  Fig.  2(a);  (2)  A  kidney  phantom  region 
registered  to  simulated  MRI,  shown  in  Fig.  2(b) 
(multimodal  registration);  and  (3)  The  same  as  (2),  but  using 
biorthogonal  wavelet  filters  to  smooth  the  B-scan  region. 

In  2D-2D  registration  experiments,  the  low  SNR  kidney 
phantom  region  was  registered  to  the  simulated  MRI  image. 
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Fig.  3.  (a)  3DUS  heart  volume,  (b)  Extracted  slice. 


2D-3D  registration  was  performed  on  US  heart  images. 
From  a  148xl60xl41-voxel  volume,  shown  in  Fig.  3(a),  a 
64x64-pixel  slice  (Fig.  3(b))  was  extracted  at  an  arbitrary 
orientation  and  registered  with  the  volume.  This  provides  a 
“gold  standard”,  as  the  correct  registration  is  known,  and 
exact  values  are  sought.  Five  trials  were  performed  with  10 
initial  points  in  3D  for  2D-2D  registration  (translation  along 
the  x  and  y-axes,  rotation  about  the  ’-axis)  and  6D  for  2D- 
3D  registration  (translation  along  the  x,  y,  and  z-axes, 
rotation  about  the  x,  y,  and  7-axes),  at  varying  distances 
from  the  true  point.  The  optimization  algorithms  were  PDS, 
the  tabu  hybrid,  and  adaptive  simulated  annealing  with  PDS. 

For  all  experiments,  a  registration  is  classified  as 
“correct”  if  the  Euclidean  distance  of  the  final  translation 
from  the  true  reference  translation  is  less  than  0.25  pixels 
and  if  the  Euclidean  distance  of  the  rotation  angles  is  less 
than  0.25°.  Optimization  was  performed  as  minimization  of 
the  negative  of  the  similarity  metric. 


VI.  Results  and  Discussion 

The  MI  and  CR  surfaces  are  shown  in  Figs.  4(a)  -  4(h). 
MI  is  characterized  by  sharp  narrow  optima.  The  CR 
functions  are  smoother,  but  also  contain  many  local 
extrema.  Filtering  the  B-scan  slightly  increased  the  CR 
maximum,  but  had  little  effect  on  the  maximum  MI  value. 

Percentages  of  correct  registrations  for  2D-2D  and  2D- 
3D  are  shown  respectively  in  Tables  1  and  2.  Registration 
accuracy  clearly  decreases  with  increasing  distance  of  the 
initial  point  from  the  true  point.  Registration  performed 
using  only  PDS  on  MI  was  very  accurate  for  initial  points 
close  to  the  reference,  but  less  accurate  for  the  CR  metric. 
This  result  is  expected,  as  CR  functions  have  more  local 
extrema  into  which  direct  optimization  can  get  trapped. 

The  tabu  hybrid  had  a  higher  percentage  of  correct 
registrations  for  all  experiments  than  PDS  when  the  initial 
point  is  far  from  the  reference,  although  these  values  are  still 
relatively  low.  For  2D-3D  registration,  the  tabu  hybrid  with 
CR  produced  the  best  overall  results,  although  PDS  was 
better  for  close  initial  points.  ASA  with  CR  did  not  have  a 
good  average  accuracy,  but  it  performed  better  than  PDS 
and  the  tabu  hybrid  for  distant  initial  points.  ASA  with  MI 
was  the  worse  performer.  The  narrow  optimum  in  the  MI 
functional  is  a  plausible  explanation.  However,  the  tabu 
hybrid  with  MI  was  much  better  than  ASA,  indicating  that 
the  former  more  thoroughly  samples  the  search  space. 


y-trans. £ 


_ 

* 


(g)  (h) 

Fig.  4.  Similarity  metrics  for:  (a-b)  MI/CR  for  low  SNR  phantom  and 
simulated  MRI.  (c-d)  MI/CR  for  low  SNR  phantom  smoothed  with  a 
wavelet  filter,  (e-f)  MI/CR  for  phantom-to-phantom  registration  at  0° 
rotation,  (g-h)  MI/CR  for  phantom-to-phantom  registration  at  12°  rotation. 


It  is  observed  that  MI  outperformed  CR  for  2D-2D 
registration,  but  CR  was  the  top  performer  for  2D-3D 
experiments.  This  result  is  expected,  as  the  heart  volume 
has  a  relatively  high  SNR  with  low  speckle,  so  that  the 
functionals  are  consistent  with  those  in  Figs.  1(a)  and  1(b). 
For  CR,  there  is  a  prominent  and  broad  optimum  with  few 
local  extrema,  while  the  MI  optimum  is  extremely  narrow. 
In  the  multimodal  2D-2D  experiments,  local  optima  are 
present  in  both  functionals,  but  those  of  CR  are  smoother 
and  broader,  and  can  be  mistaken  for  promising  regions. 

The  number  of  cost  function  (MI/CR)  evaluations  was 
also  analyzed,  and  the  mean  values  are  shown  in  Table  3. 
As  PDS  examines  only  one  region,  it  gives  the  fewest 
evaluations.  Overall,  successful  registrations  with  the  tabu 
hybrid  method  required  fewer  evaluations  than  ASA. 
Although  the  numbers  of  evaluations  are  relatively  high, 
parallelization  can  greatly  improve  the  efficiency  of  the  tabu 
hybrid,  which  is  inherently  parallel. 
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Table i 

CORRECT  2D-2D  REGISTRATIONS  FOR  OPTIMIZATION 
TECHNIQUES.  DISTANCES  ARE  MEASURED  IN  PIXELS, 
_ ROTATIONS  IN  RADIANS  (5  TRIALS). _ 


Dist. 

Dist. 

Percentage  of  Correct  Registrations  | 

(trans) 

(rot) 

PDS 

PDS 

TH 

TH 

ASA 

ASA 

(MI) 

(CR) 

(MI) 

(CR) 

(MI) 

(CR) 

4.47 

0.30 

100 

100 

100 

100 

40 

40 

7.00 

0.00 

100 

100 

100 

100 

20 

40 

13.00 

0.20 

100 

0 

100 

100 

40 

80 

14.14 

1.10 

100 

0 

80 

100 

60 

20 

14.87 

0.20 

100 

0 

60 

40 

40 

60 

15.13 

0.40 

0 

0 

60 

60 

0 

60 

17.00 

2.10 

0 

0 

40 

20 

40 

40 

35.34 

0.40 

0 

0 

60 

0 

40 

40 

42.04 

0.20 

0 

0 

40 

20 

40 

80 

61.85 

2.10 

0 

0 

0 

0 

60 

60 

Table  n 

CORRECT  2D-3D  REGISTRATIONS  FOR  OPTIMIZATION 
TECHNIQUES.  TRANSLATION  DISTANCES  ARE  MEASURED  IN 
_ PIXELS,  ROTATIONS  IN  RADIANS  (5  TRIALS). _ 


Dist. 

(trans) 

Dist. 

(rot) 

Percentage  of  Correct  Registrations 

PDS 

(MI) 

PDS 

(CR) 

TH 

(MI) 

TH 

(CR) 

ASA 

(MI) 

ASA 

(CR) 

1.41 

0.17 

100 

100 

100 

100 

0 

40 

3.32 

0.24 

0 

100 

80 

100 

0 

60 

3.46 

0.07 

100 

100 

100 

100 

0 

40 

4.12 

0.19 

100 

100 

60 

80 

0 

20 

5.83 

0.21 

0 

100 

60 

80 

0 

20 

8.83 

1.18 

0 

0 

0 

80 

0 

40 

11.79 

0.39 

0 

0 

60 

20 

0 

20 

18.25 

0.29 

0 

0 

20 

80 

20 

40 

29.88 

0.40 

0 

0 

40 

0 

0 

40 

60.61 

1.75 

0 

0 

0 

20 

0 

20 

Table m 

MEAN  NUMBER  OF  COST  FUNCTION  EVALUATIONS  (xlOOO)  FOR 
SUCCESSFUL  REGISTRATIONS. 


Type  of 

PDS 

PDS 

TH 

TH 

ASA 

ASA 

Reg. 

(MI) 

(CR) 

(MI) 

(CR) 

(MI) 

(CR) 

2D-2D 

0.316 

0.472 

1.296 

1.075 

3.261 

3.362 

2D-3D 

0.777 

0.724 

2.405 

2.606 

6.759 

7.087 

VII.  Conclusion 

This  study  demonstrates  that  US  registration  accuracy  is 
sensitive  to  similarity  metric,  optimization  method,  and 
initial  point.  The  experiments  show  that  the  tabu  hybrid 
optimizer  can  potentially  be  used  in  intra-  and  multimodal 
US  registration,  and  that  stochastic  methods,  including 
ASA,  combined  with  direct  approaches,  can  improve 
registration  accuracy.  User  interaction  in  selecting 
acceptable  starting  points  is  also  shown  to  be  important. 
Performance  is  expected  to  improve  when  these  methods  are 
integrated  with  multiresolution  techniques,  with  feature 
extraction,  or  when  used  with  combined  spatial  or  gradient 
information  [5].  The  inherent  parallellism  of  the  tabu 
hybrid  can  be  exploited  to  run  on  multicomputer  systems. 

Although  existing  direct  and  stochastic  methods  are 
effective  in  many  situations,  the  continuous  tabu  search 


hybrid,  along  with  other  stochastic  and  hybrid  techniques, 
can  be  used  as  a  robust  optimizer  for  ultrasound  image 
registration,  and  further  work  on  these  techniques  is 
warranted,  especially  for  multimodal  registration. 
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